Stochastic calculus: 
application to dynamic bifurcations 
and threshold crossings 

; 

0\ '■ Kalvis M. Jansons 1 and G.D. Lythe 2 ' 3 

To appear in Journal of Statistical Physics 

^5 \ department of Mathematics, University College London, Gower Street, 

London WC1E 6BT, England 
2 0ptique nonlineaire theorique, Universite Libre de Bruxelles CP 231, 
O ■ Bruxelles 1050, Belgium. 

o ■ 
r» : 
o . 
■ 

ctn ; 
o ; 

& ' Abstract 

For the dynamic pitchfork bifurcation in the presence of white noise, 
, ^ ' the statistics of the last time at zero are calculated as a function of the 

noise level, e, and the rate of change of the parameter, /i. The thresh- 
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cortical neuron is considered, concentrating on quantities that may be 
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1 Introduction 



Differential equations have long been used to model the dynamics of physical 
systems. With the availability of computers, the tendency to focus only on 
analysis of linear equations is being replaced by a methodology that profits 
from a judicious mixture of numerical generation of paths, bifurcation theory 
and asymptotic analysis. However, when random perturbations (i.e. "noise") 
play an important role, this new spirit is not so widespread. One reason is 
that the mathematical tools appropriate for describing stochastic paths are 
not sufficiently well known. 

In Section 2 we introduce our notation, concentrating on the aspects of 
stochastic calculus [1-6] that may be unfamiliar to applied mathematicians. 
In Section 3 we consider the dynamic pitchfork bifurcation, which is of con- 
siderable interest in its own right [8-12], and serves as an example of a system 
where noise of small magnitude has a disproportionate and simplifying effect 
[|i~6| , |i"8| . In Section 4 we discuss some exact results for threshold crossing 
problems that may be accessible to experiments but have so far received lit- 
tle attention. Our motivation comes from simple models for a single neuron's 
membrane potential pi"| , |25j and our intention is to introduce new ways of 
extracting information from numerical or experimental data. 

A strong point of the stochastic approach is the close relationship be- 
tween numerics and analysis. Throughout, we compare our calculations with 
numerical results, obtained using the algorithm described in Appendix 1. 



1.1 Dynamic bifurcations 

When used to model physical systems, differential equations have parameters 
representing external conditions or controllable inputs. Bifurcation diagrams 
are used to depict the ultimate behaviour of their solutions as a function of 
these parameters; a critical value is a point on a boundary in the space of 
parameters separating areas with different ultimate behaviour. It is natural 
to speak of a bifurcation "taking place" at a critical value. Unfortunately, it is 
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also common to speak of the qualitative behaviour of the system "changing as 
a parameter passes through a critical value" , although the paths of the non- 
autonomous system of differential equations obtained by letting a parameter 
become a function of time can differ markedly from those of the equations 
with fixed parameters [P, j9[ [I0[|. 

Performing the analysis of Section 3, we have in mind the following situa- 
tion. Looking for a symmetry-breaking bifurcation, an experimentalist slowly 
changes a parameter until a qualitative change is seen. The change is seen 
after the parameter value at which, mathematically speaking, the bifurcation 
takes place; the bifurcation is said to be delayed || |9|, |i~0l . Small random 
disturbances need to be considered explicitly to evaluate the magnitude of 
the delay. Previous work [13-17] on the dynamic pitchfork bifurcation has 
centered on the question: when do paths exceed a certain threshold distance 
from the former attractor? We calculate instead the last time at which this 
distance is zero. The cases of small and large noise are treated, and an ap- 
proximate formula is presented, with the appropriate limiting behaviors, that 
also covers the intermediate region. 



1.2 Threshold crossing and neuron dynamics 

Animal brains are made up of very many neurons, that communicate by 
firing electrical signals [24-29]. A neuron 'fires' when its membrane poten- 
tial exceeds a threshold. This membrane potential, the voltage difference 
between the interior of the cell and its surroundings, can be measured as a 
function of time. It is driven away from its rest state by many thousands of 
excitatory and inhibitory stimuli from other neurons and is sensibly modeled 
by a stochastic process. We have in mind a situation where the SDE govern- 
ing the neuron is not known in advance |25 |. We examine the dependence 
of experimentally measurable quantities on the form of the SDE describing 
the paths until they reach a threshold, concentrating on concepts that are a 
natural part of stochastic calculus but not of traditional methods. The hope 
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is that measurement of these quantities will enable the underlying dynamics 
to be deduced and lead to a greater understanding of single neurons and 
neural networks. 

In Section 4, we consider the threshold crossing problem for a fairly gen- 
eral (non-linear) scalar SDE, but with constant threshold and noise inten- 
sity. Those sufficiently well-acquainted with the probability literature will 
not find new theoretical results; our purpose is to present exact results that 
are amenable to experimental and numerical verification, despite depend- 
ing on subtle properties of the paths on small scales, and to illustrate by 
construction that exact results yield tractable expressions even when the un- 
derlying SDE is nonlinear. In Section 4.1 we split paths up into a series 
of excursions. This provides, for example, a sensitive test for asymmetries 
in the drift. In Section 4.2 we consider the amount of time paths spend in 
intervals before reaching the threshold. The function obtained in the limit 
where the intervals shrink to points is the local time; for a stochastic path, 
this is a continuous function. Section 4.3 deals with the conditioning that 
comes from considering only one part of a path and includes the calculation 
of the density of the last time at zero. 
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2 Stochastic calculus 



A stochastic differential equation (SDE) is written in the following form: 

dX t = /(X t ,t)dt + e(X t ,t)dW t . (1) 

Here X is a real-valued stochastic process; its value at time t is a real- 
valued random variable denoted by X t . If e(x, t) is always zero then ([[J) 
is just an ordinary differential equation; if e(x,t) is a real-valued function 
independent of x then (HD is a differential equation with additive white noise. 
The Wiener process W, also called standard Brownian motion, satisfies ([I]) 
with f(x,t) = and e(x,t) = 1. Its paths are continuous with probability 1. 
For any t, At > 0, the random variable W t+ At — Wt is Gaussian with mean 
and variance At, and successive increments are independent. 

Random variables are in bold type in this work. For any real-valued 
random variable a, the probability that a < x is an ordinary function of x 
denoted by V[sl < x\. The density of a (if it exists) is the derivative of this 
function with respect to x: 

R a (x) = -^P[a < x). (2) 

This is sometimes put as: R a (x)dx is the probability that a lies in (x,x + dx). 
We maintain the notation (0), with the random variable a as a subscript, 
except for the density of X t , when we write simply R(x). 

Solving an SDE on a computer consists of generating a set of values: 
{Jf (ij); i — 1, . . . , n}, where < t\ < ti < . . . < t n — t, that approximate the 
values taken by one path of X at the times tj. In the lowest order algorithm 
X(ti + i) is generated from X(ti) by adding a deterministic increment and a 
random one: 

X(t i+1 )=X(t i )+f(X(t i ),t i )At + e(X(t i ),t i )^At i n i , (3) 

where At, = t i+ i — ti and each is independently generated from a Gaussian 
density with mean zero and variance 1. Properties of the ensemble of paths, 
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such as (X t ), the mean value of X 4 , are estimated by repeating this procedure 
as many times as necessary. 

In the SDE ([!]), the coefficient f{x,t) is the mean displacement or 'drift': 
if X( = x then 

f(x,t) = ton — (X t+ At -x). (4) 

The second term on the right-hand-side of (p]) and (^) does not enter (f|) 
because it has mean zero. However, if we consider the mean squared dis- 
placement then 

e 2 (x,t) = |im i-(( Xt+At -a;) 2 ). (5) 

Further insight is gained by making the following construction. Define sets 
of times {tf, % = 0,1, ... ,1} such that = t < t\ < . . . < ti — t and 
t i+ i — ti — t/l. Choose any path of X and let 

[X] t = lim £(X, i+1 - X ti ) 2 . (6) 

1 = 

Then [X] , called the quadratic variation of X, satisfies 

d[X] t = e 2 (X t ,t)dt. (7) 

(This is true of any path of X with probability 1.) Note that if e is indepen- 
dent of X t then (|?p is an ordinary differential equation. The informal way to 
understand quadratic variation is to write 

(dW t ) 2 = dt. (8) 

Stochastic calculus [1-6] is the calculus of (continuous but non-differentiable) 
paths for which the limit (^) is non-zero. 

If X and Y are stochastic processes obeying SDEs of the form (|I|), then 
the Ito integral, 

I t = / 4 Y s dX s , (9) 
is itself a stochastic process given by the following limit: 

l t = jim^Y ti (X ti+1 -X ti ). (10) 

i=0 
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Its quadratic variation is 

pq t = rY s 2 d[x] s . (ii) 

Jo 

As suggested by (§), second order infinitesimals are not always negligible 
in stochastic calculus. This is reflected in the Ito formula, which is the chain 
rule of stochastic calculus. The SDE for /(X), where / is a C 2 function, is 
related to that for X by [1-6] 

d/(X t ) = /'(X t )dX, + A/"(X t )d[X] t . (12) 

An alternative definition of the stochastic integral is sometimes conve- 
nient, in which Y t . in ([TDD is replaced by Y tg , where t s = |(tj + t«+i). If 
this alternative definition, called the Stratonovich integral, is chosen then 
changes of variable can be performed without the extra term that appears 
in (|T^). (An extra term appears instead in the numerical algorithm; in the 
Stratonovich convention, @ no longer corresponds to the SDE ([!]).) 
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3 The dynamic pitchfork bifurcation 



One of the most common transitions observed in nature is from a symmetric 
state to one of a set of states of lower symmetry jr2| . Here we concentrate 



on the example of a pitchfork bifurcation, described by the normal form pi 

x — gx — x 3 . (13) 

If the parameter g is less than the solution x = is stable; if g > it is 
unstable and the two states x = ^fg and x = —y/g are stable. That is, for 
g > 0, any small perturbation away from x = will result in the system 
breaking the reflection symmetry by selecting one of these states. In this 
work we examine what happens if g is slowly increased through with white 
noise as the perturbation. 

In a real experiment or in a numerical search of the parameter space of 
a set of differential equations, g will typically be a complicated function of 
controllable parameters. Imagine that the experiment begins with g < and 
that the inputs are changed in such a way that it slowly increases. In the 
normal form (|13"D, x has already been scaled so that the coefficient of the 



cubic term is unity; we also rescale time so that the starting value of g is — 1. 
We thus solve the SDE 

dX t = (fitX t - X t 3 )dt + edWt, e, /i > 0, t > t = — , (14) 

fi 

and derive the statistics of the last time, s, at which X t = 0. For convenience, 
we take initial condition X to = 0. However, we shall see that the presence of 
noise removes the dependence on initial conditions for /i| loge| < 1. 

Two paths of flTil) are shown in Figure [I], where g = fit. The first path 
lingers near X 4 = until g > - /p, jumps away rather suddenly and does not 
return to 0. Although the path looks smooth on the scale of the figure, the 
moment of this jump is in fact controlled by the magnitude of the noise. The 
characteristic value of g for the jump [13-16], 



g = V2//|loge|, (15) 
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Figure 1: Dynamic pitchfork bifurcation. Two paths of the SDE ( 14|) are 
shown. In (a) noise is only important when X t is near zero. The time of 
the jump away from zero is calculated from the solution of the linearized 
equation. Once one branch of the pitchfork is chosen, a further crossing of 
is improbable. In (b) we treat the problem as one of jumping between 
minima of a potential that is slowly deepening. 
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is also found in the Hopf and transcritical bifurcations, and in the analogous 
stochastic partial differential equations with a time-dependent critical pa- 
rameter |20] . These bifurcations share the property that a fixed point ceases 



to be an attractor, but continues to exist, after the bifurcation [22]. On the 



other hand, the dynamic saddle-node bifurcation, where a fixed point ceases 



to exist at the critical point, has a characteristic delay that scales as /i» [23 



The result ( |To] ) is obtained by estimating the value of g at which X t is 



large enough that the linearized version of (|T4]) is no longer valid. In Section 
3.1, we calculate the distribution of s using the same method and show that 
it produces a good approximation for e -C ■JJl. We also give an expression for 
the value of g at which the distribution of X 4 becomes bimodal. In Section 3.2 
we consider the situation, illustrated in Figure 1(b), where paths move back 
and forth between branches of the pitchfork numerous times before settling 
down to one branch. Here we derive the density of s as an integral over the 
standard approximation for the rate of crossing between the two branches. 
This yields a good approximation for the density of s when e ^> JJl. We 
close Section 3 with a formula that interpolates between the two limits. 

3.1 Linearized equation: small-noise case 

The solution of the linearized version of (|H|) (without the cubic term) is 

X t = e^* 3 "® (x t0 + e jf* e-2^ 2 -*o)dW s ) . (16) 
That is, X 4 is a Gaussian random variable with mean 

(X t )=eI^ 2 -*o)X to , (17) 

and variance 

(X t 2 ) - (X t > 2 = v(t ,t) = eV* 2 fe-^ds. (18) 

J to 



Suppose that X t = x at some time t. Let t/, be the first time after t at 
which X is 0, with = oo if there is no such time. Note that, according to 
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fll^l), X t — > ±00 as £ — > 00. Because of the reflection symmetry, if the path 
is at zero at some time then it has an equal probability of being positive or 
negative at a later time. Thus the probability of never crossing zero after 
time t is given by 

^[Xoo > 0\X t = x} = \V\t h < oo|X t = x]+ V[t h = oo|X t = x\. (19) 

Since V[t h = oo\X t = x] + V[t h < oo|X t = x] = 1, we find 

V[t h = oo\X t = x]= 2V[X oa > 0|X t = x] - 1. (20) 

We now evaluate the right-hand-side of (B^) using (|T7|) and (ET^) . If X 4 = x 
then, for any T > t, Xt is a Gaussian random variable with mean m(T) = 
e 2 M ( T ~* and variance v(t,T). Thus 

/•m(T) 1 / 2 \ 

2P[X T >0|X t = x]-l = 2y o (27rt;(£,T))-2expf-^-^jd Z /. (21) 



Now making a change of variable in the exponent of (|2T| ) and taking the limit 
T — > 00 gives 

P[t fc = oo|X t = ar] = erf f^e"^ (2 e^ds)" 2 j (22) 

where erf (it) = ^= /J 1 e~ y2 dy. 

The last time at zero is the random variable s. To calculate its density, 
we replace the initial condition x by a Gaussian random variable with mean 
and variance given by (|18|). Then Hi(t) = V[t > s] is obtained by integrating 
([22] ) over the Gaussian density: 

/■oo 1 

HAt) = 2/ P[t A = oo|X t = x](27rv(t))~2 exp 







/ e"" erf 


!■( 


/o 





7 t °°exp(-/is 2 )d g y 

j/ o exp(-/XS 2 )ds / 




(23) 



^ I V /t exp(-/xs 2 )ds / 



11 



The density of the last time, s, at which X t passes through is 

Ut) = Ai/ l(t ) = I- ^t^l p (24) 

(/to exp(-/is 2 )ds J t °° exp(-/is 2 )ds) 2 



For fi — > 0, this density is symmetric about t = with width proportional 
to -4-. The above calculations do not include the effect of the cubic term in 
Because this term always pushes paths towards X t = 0, the function 



H\{t) from fl23|) is an upper bound on the probability of no crossing of zero 
after time t. 

Self-consistency of the approximation (|2~3l) demands that X t 3 <C gX 4 for 
t < From ([18]), X t ~ e for t < Thus we require e <C ^//L The 
calculations of this Section are quantitatively accurate when, in addition, 
the probability of a crossing of zero for t ^> -k= is negligible, corresponding 
to the situation shown in Figure 1(a). This assumption will be examined in 
the next subsection; it is also true for e <^ ^fji. 

We have used a zero initial condition in (|i~4l) because it is convenient 
to have (X t ) = 0. This initial condition corresponds to the case where 
the system is already close to its equilibrium before the sweep is started. 
With a non-zero initial condition, our calculations remain valid if (X t 2 ) > 
(X t ) 2 for t > 0. Comparison of the magnitudes of ( |17D (proportional to 
X to ) and flT8| ) (proportional to e) shows that a non-zero initial condition is 
unimportant unless X 4o > eexp(^j). This insensitivity to initial condition 
is characteristic of the simplification produced by noise in problems where 
paths sweep repeatedly past an invariant manifold |L6], [LB]. 

The results of this Section are changed little if the noise is 'colored' (i.e. 
if W t in (|T4]) is replaced by a stochastic process whose increments are corre- 
lated) | pLQ|1 . If a deterministic perturbation is used instead of noise, the effect 



is proportional to the mean of the perturbation. Multiplicative noise is less 
important than additive noise |L9]. gOj . 
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3.1.1 Transition to bimodal density 

To estimate when nonlinear terms become important in a dynamic bifurca- 
tion, one can calculate from the linearized version the time when X t exceeds 
some threshold [13-16]. Here, we calculate instead, from the nonlinear equa- 
tion with e <C y/JI <C 1, the time when the distribution of X 4 changes from 
unimodal (one maximum at 0) to bimodal. We do this by considering the 
evolution of an ensemble of paths from a Gaussian distribution of initial 
conditions into the latter part of the path where |X t | — > ^fg from below. 
Effectively, noise provides a random initial condition for the subsequent, de- 
terministic, evolution. 

The solution of the deterministic equation, 

x = fitx — x 3 , (25) 

with initial condition x = Xi at t = t\ is: 

x = ^-f^~ tl) (V + 2 [ e^^ds) 5 . (26) 

Now, instead of a constant as initial condition, we take the following Gaussian 
random variable: 

/ 7T \ 3 I t 2 

X tl = w(ti)n, where w{t) = e I - J e2 M (27) 

is the large-t limit of v(t) flUf), and n is a Gaussian random variable with 
mean and variance 1. We use 2 f t exp(/xs 2 )ds ~ exp( / ut 2 )( t ut)~ 1 and obtain 

^y^) 1 - <28) 

The density of X t is then given, with t as a parameter, by 

R( X) = f- w(t ) exp (_^ pM ) ( X " S) ' ' <29) 
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it has a maximum at a non-zero value of x if w 2 (t) > The time when 
the density changes to bimodal therefore satisfies 



Qt)* - ^log if = 2/x| loge| + | log 1. (30) 
For small /i, we recover the result {lit) 2 = 2/i|loge| that comes from the 



analysis of the linearized system with a fixed threshold |13|, [14], [15], |16 |. 



3.2 Effect of the cubic term: large-noise case 

The probability that there is no crossing of zero after time t is a positive 
increasing function of time, H(t). We can therefore make the decomposition 

±H(t) = r(t)H(t). (31) 

So 

H(t) = exp (- J™ r(s)ds\ . (32) 

The function r(t) is everywhere positive and can be interpreted as the proba- 
bility per unit time of a crossing of zero. In Section 3.1 we obtained an upper 
bound on H(t) from the linearized equation. A lower bound on H(t) can be 
obtained as follows. Calculate the probability per unit time of a crossing 
of zero from flli]) with fit taken as a fixed parameter. Denote this function 
by r 2{t)- Because the slow increase of g makes crossings of zero rarer as t 
increases, r 2 (t) is an upper bound on the rate of zero crossings. The function 
H 2 (t), calculated from (^), is therefore a lower bound on the exact result 
H(t). For direct comparison with numerics, Hzif) is obtained by integrating 
T2{t) over the density of X t ; the result is an unwieldy expression. Here we 
restrict ourselves to the case g > e and examine the validity of the compact 
expression that results. 

The expression ( |3T]) is most useful and intuitive for g > e, when paths 
spend most of their time near ±y/g. Then r 2 (t) is well approximated by the 
inverse of the mean passage time from X t = to X t = [1], [22]] : 

rtf) = ^g exp (-Q . (33) 
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Using (|3"3"| ) in (RJ) with g = fit gives 



2 



*M = ^«p(-^§«p(-£JJ. < 34 > 

The expression for the density of s corresponding to (|34]) is 

-ex (-^ 



/? s (/) = Aff 2 (f) = ^ exp (_£.) exp f-^-exp f-|l ) } . (35) 



As /i — > 0, the mean value of s is |L6 



2\ \ 2 



<s) = K s * + v) where 9 " = £ ( 21oe (^tJJ (36) 

and Euler's constant 7 = 0.57721 .... In the same limit, the standard devia- 
tion is 

The approximations (|33| ) -(|3~7| ) are valid if g* ^> e, which is true if e > ^/Ji 



^)-^) 2 Y = ^=A^ g ^-\ . (37) 



Note also that H 2 (t) ~ 1 for (7 > ^//I if e < y'/Z, as required for validity of 
the calculations of Section 3.1. 

The decomposition (|3l| ) also provides a natural way to unify the above 



approximation with that of Section 3.1. Let r(t) = r\(t) +7*2 (t) where r\(t) 



Hi(t) 1 4:Hi(t), with Hi(t) given by (^). Then H(t) can be obtained from 



([32] ) . This yields the following approximation to R s (t) that has ( p5[) and ( 23 ) 
in the large-noise and small-noise regimes: 

R s (t) = (r 1 + r 2 )H 1 H 2 

= H^H^+H^H^t). (38) 

In Figure ||, we used 

fsexp(-| 
^eexp(-i) 



,wH?!!}t j (39) 
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Figure 2: Last crossing of zero: numerical results vs formula. The dots 



are the histogram compiled from numerical generation of 10000 paths of fll4|) 
and the solid lines are the uniform approximation described in Section 3.3. 
Figure (a) exhibits the small- noise case (e < y 1 'ji) and (b) the large- noise 
case. Note g = fis. 
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4 Exact results on threshold crossings 



In this Section we consider processes obeying SDEs of the following type: 



dX t = fpQdt + dW t 



(40) 



with initial condition Xo = 0. Paths are followed until the first time, t&, 
that they reach a threshold b > 0. If f(x) is not known in advance, what 
can be learned from the paths of X t for t < t&? The process X t models, 



for example, the membrane potential of a neuron |24], [30] or the path of a 
diffusing particle |3Tfl . To produce (pUf ) from the more general SDE ([I]), we 
assume that e is constant. We can then scale time so that the coefficient of 
the second term in (|1|) is 1. We also primarily consider the case where the 
drift, f(x), is towards x = 0. 

We do not calculate two of the quantities traditionally considered by theo- 
rists: the density of X t and that of t b p7|, ^8f . We calculate instead quantities 
that are measurable by an experimentalist who can record individual paths, 
but have so far received little attention outside of probability theory All 
quantities calculated are derived from the following function: 



S{x) 



f X exp{2U(y))dy x > 0, 
Jo 

[°exp(2U(y))dy x<0, 

i J X 



(41) 



where 



U{x) 



f{y)dy- 



(42) 



The function S(x) is generated on a computer by solving the ordinary differ- 
ential equation 



S'(x) = exp(2U(x)) x>0, 
5(0) = 0, 

S'(x) = -exp(2U(x)) x<0. 



(43) 
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Further, S(x) has the following simple interpretation. Given X t = x, < 
x < b, at some time, the probability that X reaches b before is 
An advantage for performing calculations is that S(x) is time-independent. 
This is a consequence of the 'memoryless' property of paths of the SDE: the 
future of X t depends only on the value of X 4 , not on the path before this 
time. In the traditional approach, based on the Fokker-Planck equation, 
measurable quantities are derived from the density of X t , which is time- 
dependent because of the flux of probability through the threshold. 

A path of a stochastic process obeying the SDE (^) can be thought of 
as being made up of a series of excursions from zero. Each time X t returns 
to the process is considered as starting afresh, and each excursion treated 
as an independent event [2-5]. In Section 4.1 we derive the statistics of the 
excursions in terms of the maximum distance from the origin during the 
excursion. 

In Section 4.2 we use the concept of local time, which measures the oc- 
cupation density of a path [|, |, |, [3^| . Given a path of a stochastic process 
satisfying (f40|), the local time at a is a stochastic process whose value at time 
t is given by 

1 r* 

L t (a) = lim — — / 1{X S G (a — Aa, a + Aa)}ds, (44) 
Aa^o 2Aa jo 

where X{0} = 1 if is true and is zero otherwise. Note that for fixed a, 
L 4 (a) increases at times when X t = a and is constant otherwise. 

Local time is the limit of functions that can be constructed from a series 
of values of X t that could be available to an experimentalist. Suppose the 
amount of time spent in boxes (a — Aa, a + Aa) before tb is recorded. The 
function obtained in the limit Aa — ^ is L ti) (a). In Section 4.2 we consider 
two ways of comparing numerical or theoretical data with L tfj (a): as an 
average over many paths and for one path at a time. 

In Section 4.3 we consider the portion of the path of X after s&, the last 
time that the path passes through zero before tj,. This is convenient, both 
numerically and experimentally, when t;, is large: rather than generating 
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Figure 3: Threshold crossing. The first time a path attains a level b is a 
random variable, denoted by t b . Also indicated is s b , the last time before t b 
at which the path is at 0. The path shown is a path of the Wiener process. 

or recording the whole path, one can concentrate on only the last portion. 
Further, it is one case where calculations based on a constant threshold will be 
a good approximation to the case where the threshold approaches a constant 
value for large t. 
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4.1 Excursions 



An excursion is a portion of a path between successive crossings of 0. A 
path can be thought of as consisting of a series of excursions occurring one 
after the other, with properties chosen according to a certain probability. We 
construct this probability as follows. For x > 0, let N(x) be the number of 
excursions before s& (see Figure ^) during which the maximum of X is greater 
than x; for x < 0, let N(x) be the corresponding number of excursions with 
minimum less than x. 

Suppose that at some time t, X 4 = x > 0. Then the probability that X 



reaches b before is |M. Thus 



x>0. (45) 



Now consider x < and suppose that at some time t, X t = 0. Then the 
probability that X reaches x before b is gT^gm ■ Thus 

"M*> ^ n| ^ s(#nw WW " x<0 - (46) 

From (|4*5D and (M) we find the following exact result for the mean of N(x): 

x < 0, 

(X(.r)) = { x x = (47) 
- 1 x > 0. 




We find good agreement between the exact result (47) and averages over 
paths of generated with a finite time-step (see Figure f|). Note that, if the 
drift f(x) is symmetric, then (N(—b)) = 1. This is a sensitive test for 
asymmetries in f(x). We have not counted the portion of the path after Sf, 
as an excursion. We return to this in Section 4.3. 
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Figure 4: Pre-threshold excursions. Here (N(x)) is the average number of 
excursions with maximum greater than \x\ before t& (6 = 2). The curves, 
produced by averaging over 100 numerically-produced paths with At = 10~ 4 , 
are in excellent agreement with the exact result ([47]). 
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4.2 Occupation density 

Another way to extract information from a path followed until it reaches a 
threshold is to construct the occupation density: the amount of time spent in 
boxes (a — Aa, a + Aa) before t&. The function obtained in the limit Aa — ► 
is Lt b (a), the local time for a path run until it first reaches the threshold b. 
For convenience, in this Section we denote Lt 6 (a) simply by L Q . 

We first consider the local time, averaged over many paths. Given X = x, 
c < x < b, let l(x) be the mean value of the local time at a before the first 
time that X t exits the interval (c, b). Using the delta-function notation to 
rewrite (f44|), 

L a = / tb 5(X s -a)d[X] s , (48) 
J o 

we find (see Appendix 2) that l(x) satisfies: 

- 5(a) = f{x)l'{x) + \l"{x), 1(c) = 1(b) = 0. (49) 

With x = and c = — oo we obtain (L a ) from the solution of (^9|): 

( 2exp(-2U(a))(S(b)-S(a)) a > 0, 
^ a/ { 2exp(-2?7(a))5(6) a < 0. 1 J 

This average is displayed for several choices of f(x) in Figure |[ 

The average (|5"0|) is an average over realizations and thus is useful when 
many paths have been recorded. Quantitative comparison can also be per- 
formed for just a few paths, exploiting the fact that an occupation density 
can be constructed from just one path. A remarkable fact is illustrated in 
Figure ^|: L a is a continuous function of a when e ^ 0. In fact, L a can be con- 
sidered as a stochastic process indexed by the space variable a. In Appendix 



3 we show that L a satisfies the non-autonomous SDE [25 



( , L 2LUw z -2(f(b-z)L z -l)dz z<b, (5i) 

2L]dW z -2f(b- z)L z dz z>b, 
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Figure 5: Average occupation density. The curves are fl50l), the analytical 
result for the average (over many paths) of the amount of time spent at a 
before the first time that a path, started at 0, reaches 6 = 2. 
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Figure 6: Occupation density. Each curve was obtained by following one 
path until it reached 6 = 2. We used Aa = 0.01 and At = 10~ 5 . Note the 
discontinuities in L a for the deterministic path. 

where z = b — a. Traces mimicking those of Figure |6| can therefore be directly 
generated by solving (|5T| ) with the initial condition Lo = 0. The drift of the 
SDE (|5l"l) ensures that ~L Z remains positive for z < b (a > 0); a path is 
followed until L z = for some z > b (a < 0). 

For a > 0, the density of the random variable L a is 




(52) 



with (L a ) given by @ 



The corresponding result for a < is 





exp( 



y 



y>0 



(53) 



(54) 



aKyj S(a) + S{b)q{a) 
where q(a) = 2 exp(-2[/(a)) (S(a) + S(b)). 



q(a) 
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4.3 Last zero crossing 

When tf, is large, it may be impractical experimentally to record the whole 
path until t&. Numerically, it also becomes time-consuming to generate 
enough paths to evaluate averages. One way around this is to concentrate 
on the portion of the path between s b , the last time that the path passes 
through zero, and t&. Numerically, it is possible to generate this portion of 
the path directly. 

Generating one portion of a path conditioned on desired properties is 
equivalent to selecting a subset of possible paths (here, paths that do not 
return to 0). This is equivalent to changing the drift of the SDE. Define 
Y to be the conditioned process and suppose X t = x at some t. Then the 
probability that t > s b is h(x) = That is, the probability that a path 
of X with X t = x is also a path of Y is h(x). Let AX = X t+A * — X t . Then 
the drift of the SDE for Y is 

r 1 /v \ v 1 (Ms + AX) AX) 

&? At {Yt+At - x) = &$>Xt (Ms + AX)) 

= hm 1 ((^) + ^)AX + ...)AX) 

At^oAt h(x) v ' 

V 1 /V \M S ' (X) 

hm — (X t+A t - x) + 



A^oAt x 1 S(x) 

Because we are choosing a subset of paths with non-zero probability, the rate 
of increase of quadratic variation (which is the same for each path) is not 
affected by the conditioning. The SDE for Y s is thus 

dY s = /(Y s )ds + dW s where f(y) = f(y) + |M. ( 56 ) 

The transformed drift f(x) is singular at x — 0, which has the required effect 
of preventing paths of Y from returning to 0. (Note the following relationship 
between / and /: 

f'( x ) = f'( x )-P( x ) + f( x ). (57) 
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Figure 7: Last zero crossing. The curves are t& — s& versus the height of 
the boundary b for various choices of the drift f(x). They are produced by 
integrating the expression 

The boundary condition for (|57|) is xf(x) — > 1 as x — > 0.) 

From the transformed drift, quantities such as U'(y) = —f(y) can be 
defined and calculated as described in this section. For example, the mean 
value of t;, — St is the mean first passage time from Y = to Y = b, which 
is given by (Appendix 2): 

(t 6 - s b ) = lim 2 J* exp (2U(yj) J* exp (-2U(z))dz dy (58) 

= lim 2 [ b -±—exp{2U(y)) F S 2 (z) exp (-2U(z))dz dy. 
x^-o J x b z (y) Jo 

In Figure ^| this mean is displayed as a function of the threshold b for several 
examples of f(x). 
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5 Conclusion 



If a symmetry-breaking bifurcation is searched for by making a parameter 
a function of time, the delay produced means that it is not obvious how to 
define or measure its position using standard ideas. The presence of noise 
simplifies the picture by removing the dependence on initial conditions in the 
limit of small sweep rate, but means that the usual criteria for determining 
the position of the bifurcation, based on deterministic differential equations, 
are no longer appropriate P2| . We study the dynamic pitchfork bifurcation 
and choose a simple concept: the last time, s, that a path is at zero. The 
density of s is calculated as a function of the noise level, e, and the rate of 
change of the parameter, \i. If e y/ji, there is a neat separation between 
the noise-dominated and nonlinear parts of the path. The mean of s in 

this case is the point of loss of stability as calculated by static theory and 

_i 

the standard deviation of s is proportional to fi 2. An expression is also 
calculated for the time when the density of X t becomes two-humped. Here 
the result is consistent with those obtained by setting a threshold on X t : the 



characteristic delay of the bifurcation is y 2/i| loge|. In the second case, e 
y/JI, the analysis is analogous to that of a two-state system. The mean value 
of s is a slowly increasing function of fi and e, and the standard deviation 
of s is proportional to e. Numerical results are compared with a uniform 
approximation that captures these two cases as limits. 

A threshold crossing problem is the ideal setting for stochastic methods. 
In Section 4, we consider some exact results and compare them with nu- 
merical data. Our aim is to calculate quantities that an experimentalist, 
confronted with a series of paths like that of Figure || can use to deduce the 
underlying dynamics. The particular motivation comes from the modeling of 
neurons that fire when their membrane potential exceeds a threshold. Results 
are expressed in terms of the function S(x) (f4~l~D, which is easily generated 
on a computer. 
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Appendices 

1. Numerical methods 

As suggested by the SDE notation (0), the simplest algorithm for solving 
an SDE is the following. At each step add two increments to X ( : a drift 
given by /(X t , t) At, where At is the timestep, and another that is a random 
variable proportional to e\fAt. The stochastic version of the Euler method 
for constructing an approximation X t to X t is: 



X t+ At = X t + f(X t , t) At + e (X t , t)VAt n (59) 

where n is a Gaussian random variable with mean zero and variance 1. The 
method that was used to generate the figures in this work is analogous to 
the second-order Runge-Kutta method 



Xt+At — X t — 

| (j(X t , t) + f (X t + f(X t , t)At + eVAl n, t + At) ) At + e^At n. (60) 

Note that only one Gaussian random variable is generated at each time step. 
For the algorithm ([SOf) there exist positive constants k and S such that, if e 
is constant and if the step size At is less than 5, then \ (h(X t )) — (/t(Xt))| < 
K(At) 2 for any polynomial function h. Higher order algorithms exist 0. 
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2. Derivation of PDEs from the Ito formula 



Consider the SDE 

dX t = /(X t ) dt + e dW t , X = x, a < x < b, (61) 

and define T(x) to be the mean value of the first time that Xf = b. Since 
T(x) is an ordinary function of the initial condition x, we can consider the 
evolution of T(X 4 ) along a path. The SDE for T(X t ) is, using the Ito formula 

(0): 

8T d 2 T 
We know, from the definition of T(x), that 



dT (X,) = — ( X *)dX t + |e 2 — (X t )d[X] t . (62) 



(dT(X t )) = -dt. (63) 

From (|7|), the quadratic variation for the process X 4 satisfies d[Xj] = e 2 dt 
and, from (||), (dX f ) = /(X t )dt. Thus taking the mean value of both sides 
of (0) gives § § 



l 2 d 2 T | &T 
dx 2 dx 



+ =~1. T(a)=T(6) = 0. (64) 



In Section 4, we also consider the quantities h(x), the probability that X 
reaches b before a, and l(x), the mean value of L tfc (a). The appropriate PDEs 
are derived as above, using (d/i(X t )) = and (d/(X t )) = — 5(a)dt instead of 
(1631). 
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3. SDE for local time in space variable 



As in Section 4.2, L a denotes the local time at a (f44|), evaluated at t = t^. We 
consider L a as a stochastic process indexed by a and derive the corresponding 
SDE. By differentiating (EOT), the drift of the SDE for L a is: 



d = 2(/( )L a - 1) x > 0, 
do N \ 2/(o)L a i<0. 

To find the quadratic variation of L a , we begin by using Ito's formula to 
derive an SDE for the function u(X t ) where X t satisfies (HO) and 



u(x)= { X X - a ' (66) 
a x < a. 

Note that u"(x) = 5(x — a). The SDE for w(X t ) is 

du(X t ) = l{X t > a}dX t + \5{X t - a)dt. (67) 
Integrating (|67D yields 

u(X t ) - u(Xo) = f X{X t > a}dX t + | f 5(X t - a)dt. (68) 
Jo Jo 

Choosing t = tb and X to = gives 

u(b) - w(0) = [ tb I{X t > a}dX t + \L tb {a). (69) 
Jo 

Thus 

L a = 2Z a + 2(m(6) - u(0)) where Z a = /*" l{X t > a}dX t . (70) 

JO 

Now define a set {<2j; i = 0, 1, . . . , /} such that a = a < a\ < . . . < ai = b 
and cii + i — ai = (b — a) /I and denote the quadratic variation of L a by [L] a . 
Then 

[L] a = hm^(L ai+1 -L ai ) 2 

= lim^4(Z 8!+1 -Z a ,) 2 . (71) 

i=0 
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To evaluate ([H]), we define the processes = ( J t I{ai < X t < a i+ i}dX^ . 
Then = (z a . +1 - Z a ^f '. The SDE for Q| is, using Ito's formula @, |2§, 

dQ? = 2 (/ Q *X{a, < X t < a i+ i}dX t ) J{a, < X t < a !+1 }dX ( 

+J{ai < Xj < a i+1 }dt. (72) 



Thus 



= 2 £ (Z ai+1 - Z ai ) X{a, < X t < a m }dX t 

+ [ b l{ ai <X t <a i+1 }dt. (73) 
Jo 

In the large-/ limit, only the contributions to [L] from the second integral 
in (|75|) remain [H, 120] and 



[L] a = 4/ tb J{X t >a}dt 



(74) 



Thus, defining z = b — a, L z satisfies the SDE 



i 



2LldW z -2(f(b-z)L z -l)dz z<b, 
dL z = ■{ i (75) 

2Ll dW 2 - 2/(6 - z)L z dz z > b, 

with initial condition Lq = and an absorbing boundary at 0. 
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